Binary effects of sediment on the mortality of adult corals
Lillian J. Tuttle
9/24/2020
Here I will explore how DEPOSITED AND SUSPENDED sediment affects the MORTALITY of tropical, scleractinean coral ADULTS. The database to be used is comprised of data extracted from studies deemed relevant during the systematic literature review process.
The Dataset
Deposited Sediment
Specifically, we will explore the results of 827 datapoints from 28 studies conducted in 3 oceans on 87 species within 46 genera.
## n
## 1 827## Ref n
## 1 DS08a 162
## 2 DS08b 113
## 3 DS71 112
## 4 DS04 70
## 5 DS68 48
## 6 DS46 44
## 7 DS29 40
## 8 DS30 36
## 9 DS17 24
## 10 DS45 24
## 11 DS16 19
## 12 DS18 16
## 13 DS27 16
## 14 DS32 16
## 15 DS07 13
## 16 DS03 12
## 17 DS43 12
## 18 DS05 10
## 19 DS33 8
## 20 DS12 7
## 21 DS26 6
## 22 DS11 4
## 23 DS21 4
## 24 DS28 4
## 25 DS25 2
## 26 DS36 2
## 27 DS69 2
## 28 DS23 1## Ref_name n
## 1 Hodgson (1989) PhD Chapter IV 275
## 2 HDR EOC and CSA Ocean Services, Inc. (2014) 112
## 3 Duckworth et al. (2017) 70
## 4 Flores et al. (2012) 48
## 5 Stafford-Smith (1993) 44
## 6 Stafford-Smith (1990) PhD Chapter 4 40
## 7 Stafford-Smith (1992) 36
## 8 Piniak (2007) 24
## 9 Rogers (1983) 24
## 10 Philipp and Fabricius (2003) 19
## 11 Piniak and Brown (2008) 16
## 12 Sofonia (2006) PhD Chapter 4 16
## 13 Stewart et al. (2006) 16
## 14 Hodgson (1990b) 13
## 15 Bessell-Browne et al. (2017a) 12
## 16 Peters and Pilson (1985) 12
## 17 Hodel (2007) 10
## 18 Vargas-Angel et al. (2006) 8
## 19 Loiola et al. (2013) 7
## 20 Sofonia (2006) PhD Chapter 3 6
## 21 Lirman et al. (2008) 4
## 22 Rogers (1979) 4
## 23 Sofonia and Anthony (2008) 4
## 24 Gil et al. (2016) 2
## 25 Shore-Maggio et al. (2018) 2
## 26 Zill et al. (2017) 2
## 27 Selim (2007) 1## Ocean n
## 1 Pacific 757
## 2 Atlantic 69
## 3 Indian 1## Selecting by n## Gsp n
## 1 Acropora_millepora 38
## 2 Leptoria_phrygia 42
## 3 Montipora_aequituberculata 26
## 4 Oxypora_glabra 23
## 5 Pavona_cactus 29
## 6 Porites_cylindrica 30
## 7 Porites_lobata/lutea 56
## 8 Porites_lutea 28
## 9 Porites_rus 30
## 10 Turbinaria_reniformis 24## Selecting by n## Updated_Genus n
## 1 Porites 152
## 2 Montipora 137
## 3 Acropora 88
## 4 Leptoria 42
## 5 Pavona 37
## 6 Turbinaria 33
## 7 Oxypora 31
## 8 Pocillopora 29
## 9 Echinopora 23
## 10 Mycedium 21Suspended Sediment
Specifically, we will explore the results of 272 datapoints from 6 studies conducted in 3 oceans on 19 species within 14 genera (shown below only the top 10 species and genera, by number of datapoints).
## n
## 1 272## Ref n
## 1 SS06 150
## 2 SS08 48
## 3 SS22c 28
## 4 SS27 16
## 5 SS22a 10
## 6 SS04 9
## 7 SS22b 8
## 8 SS25 3## Ref_name n
## 1 Browne et al. 2015 150
## 2 Flores et al. 2012 48
## 3 Rice 1984 46
## 4 Anthony et al. 2007 16
## 5 Bessell-Browne et al. 2017 9
## 6 Te 2001 Chapter 6 3## Ocean n
## 1 Indian 150
## 2 Pacific 76
## 3 Atlantic 46## Selecting by n## Gsp n
## 1 Merulina_ampliata 50
## 2 Pachyseris_speciosa 50
## 3 Platygyra_sinensis 50
## 4 Montipora_aequituberculata 30
## 5 Acropora_millepora 21
## 6 Acropora_intermedia 16
## 7 Cladocora_arbuscula 8
## 8 Manicina_areolata 8
## 9 Solenastrea_hyades 8
## 10 Siderastrea_radians 6## Selecting by n## Updated_Genus n
## 1 Merulina 50
## 2 Pachyseris 50
## 3 Platygyra 50
## 4 Acropora 37
## 5 Montipora 36
## 6 Cladocora 8
## 7 Manicina 8
## 8 Solenastrea 8
## 9 Siderastrea 6
## 10 Isophyllia 4
## 11 Phyllangia 4
## 12 Scolymia 4
## 13 Stephanocoenia 4Exploratory plots
First let’s explore all data from all species for which there is binary data about the presence of ‘small necroses’, ‘large necroses’, and ‘total mortality’ as a result of exposure to deposited sediment.
DEFINITIONS:
- ‘Small necroses’ are defined as areas of tissue mortality that encompass between 0% and 50% of the coral’s surface area.
- ‘Large necroses’ are defined as areas of tissue mortality that encompass greater than or equal to 50% of the coral’s surface area, but less than 100% (total mortality).
- ‘Total mortality’ is defined as 100% of coral’s surface area experiencing necrosis. In other words, it is the death of an entire coral entity (e.g., colony, branch, fragment, other experimental unit).
Threshold Analyses with Binary Data
Now let’s calculate the thresholds for each category of mortality, based on the binary data explored above.
DEFINITIONS:
- ‘LOAEL’, or the ‘lowest observed adverse effect level’ is defined as the lowest dose at which there was an observed adverse effect.
- ‘NOAEL’, or the ‘no observed adverse effect level’ is defined as the highest dose at which there was NOT an observed adverse effect.
- ‘Adverse effect’ is defined here as any response of a coral individual, colony, or experimental treatment group that may negatively affect the coral’s fitness and/or survival. These adverse effects may include sublethal physiological changes (e.g., significantly reduced growth or photosynthetic rates when compared to coral in ambient/control conditions), bleaching/tissue loss, and mortality. This definition of adverse effect is independent of its magnitude, so while the affect may have potentially reduced fitness, its magnitude may be sufficiently small that the fitness effect is not measureable.
DEPOSITED SEDIMENT
Small Necroses
NOAEL/LOAEL defined by exposure concentration
## LOAEL
## 1 4.9## NOAEL
## 1 4.4
NOAEL/LOAEL defined by exposure duration
## LOAEL
## 1 0.917## NOAEL
## 1 0.917
Large Necroses
NOAEL/LOAEL defined by exposure concentration
## LOAEL
## 1 20.8## NOAEL
## 1 20.8
NOAEL/LOAEL defined by exposure duration
## LOAEL
## 1 3## NOAEL
## 1 3
Total mortality
NOAEL/LOAEL defined by exposure concentration
## LOAEL
## 1 20.8## NOAEL
## 1 20.8
NOAEL/LOAEL defined by exposure duration
## LOAEL
## 1 1## NOAEL
## 1 1
SUSPENDED SEDIMENT
Small Necroses
NOAEL/LOAEL defined by exposure concentration
## LOAEL
## 1 3.22## NOAEL
## 1 3.22
NOAEL/LOAEL defined by exposure duration
## LOAEL
## 1 14## NOAEL
## 1 14
Large Necroses
NOAEL/LOAEL defined by exposure concentration
## LOAEL
## 1 29.1## NOAEL
## 1 29.1
NOAEL/LOAEL defined by exposure duration
## LOAEL
## 1 84## NOAEL
## 1 84
Total mortality
NOAEL/LOAEL defined by exposure concentration
## LOAEL
## 1 29.1## NOAEL
## 1 29.1
NOAEL/LOAEL defined by exposure duration
## LOAEL
## 1 84## NOAEL
## 1 84
Logistic Regression Model Fitting
DEPOSITED SEDIMENT
Small Necroses
## n
## 1 602## Data: smallDS2
## Models:
## modDS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS6: (1 | Gsp)
## modDS9: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS9: (1 | Ref)
## modDS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS6 5 660.76 682.77 -325.38 650.76
## modDS9 5 702.81 724.81 -346.40 692.81 0.000 0 1
## modDS12 6 652.31 678.71 -320.15 640.31 52.496 1 4.31e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: smallDS2
## Models:
## modDS10: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDS10: Ref)
## modDS11: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modDS11: (1 | Gsp) + (1 | Ref)
## modDS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS10 4 657.21 674.81 -324.61 649.21
## modDS11 5 650.97 672.97 -320.48 640.97 8.2432 1 0.00409 **
## modDS12 6 652.31 678.71 -320.15 640.31 0.6604 1 0.41641
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Gsp) + (1 | Ref)
## Data: smallDS2
##
## AIC BIC logLik deviance df.resid
## 651.0 673.0 -320.5 641.0 597
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1315 -0.6021 -0.1801 0.5730 3.0137
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 2.871 1.694
## Ref (Intercept) 2.404 1.551
## Number of obs: 602, groups: Gsp, 75; Ref, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.975421 0.562799 -3.510 0.000448 ***
## Sed_level_stand_mg 0.002889 0.001178 2.453 0.014167 *
## Sed_exposure_days 0.013832 0.004908 2.818 0.004828 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_small_necroses ~ log10_conc + log10_dur + (1 | Gsp) +
## (1 | Ref)
## Data: smallDS2
##
## AIC BIC logLik deviance df.resid
## 639.8 661.8 -314.9 629.8 597
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3280 -0.5459 -0.1809 0.5600 3.1754
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 3.219 1.794
## Ref (Intercept) 3.934 1.983
## Number of obs: 602, groups: Gsp, 75; Ref, 21
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.4867 1.1307 -4.852 1.22e-06 ***
## log10_conc 1.3053 0.3869 3.374 0.000742 ***
## log10_dur 1.5651 0.4175 3.749 0.000178 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Effect estimates
## [1] 0.8272425##
## p 0 1
## 0 344 61
## 1 43 154## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.9049## R2m R2c
## theoretical 0.07964959 0.7100653
## delta 0.07245285 0.6459074
## $Gsp
##
## $Ref
## Groups Name Std.Dev.
## Gsp (Intercept) 1.7942
## Ref (Intercept) 1.9834## [1] 7.267637
## Est LL UL
## (Intercept) 0.004141512 0.0004515197 0.03798754
## log10_conc 3.688788053 1.7279975250 7.87452360
## log10_dur 4.783110717 2.1103679520 10.84083376For every 10-fold increase in exposure concentration of deposited sediment, the odds of developing small necroses (<50% of coral tissue surface area) increase by 3.7 times (95% CI 1.7, 7.9; GLMM z = 3.374, p = 0.0007), after accounting for exposure duration and variability among studies and species. For every 10-fold increase in exposure duration of deposited sediment, the odds of developing small necroses increase by 4.8 times (95% CI 2.1, 10.8; GLMM z = 3.749, p = 0.0002), after accounting for exposure concentration and variability among studies and species.
Prediction figure
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDS11_log <- predict(modDS11_log, newdata=smallDS2, type="response")
smallDS3 <- cbind(smallDS2, pred_modDS11_log)
smallDS_plot <- ggplot(data = smallDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_small_necroses)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Small necroses due to sediment exposure?") +
scale_x_log10(limits=c(1,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDS11_log), inherit.aes=FALSE) +
theme_bw()
smallDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
By sediment exposure concentration
jvalues <- with(smallDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
smallDS2$log10_conc <- j
predict(modDS11_log, newdata = smallDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
##
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
smallDS_plot2 <- ggplot() +
geom_point(data = smallDS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_small_necroses,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of small\nnecroses due to sediment exposure",
color = "\nBinary Data by Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(1,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS1),
label=c("0.01","0.1","1","10","100","1000",round(loaelDS1,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelDS1, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
smallDS_plot2
By sediment exposure duration
jvalues <- with(smallDS2, seq(from = min(log10_dur), to = max(log10_dur), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
smallDS2$log10_dur <- j
predict(modDS11_log, newdata = smallDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_dur_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_dur","Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
smallDS_plot3 <- ggplot() +
geom_point(data = smallDS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_small_necroses,
color = Ref)) +
labs(x = "Sediment exposure duration (days)",
y = "Predicted probability of small\nnecroses to sediment exposure",
color = "\nBinary Data by Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.1,200), breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
smallDS_plot3
Large Necroses
## n
## 1 522## boundary (singular) fit: see ?isSingular## Data: largeDS2
## Models:
## modDSb6: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb6: (1 | Gsp)
## modDSb9: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb9: (1 | Ref)
## modDSb12: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb6 5 346.09 367.38 -168.04 336.09
## modDSb9 5 392.60 413.88 -191.30 382.60 0.000 0 1
## modDSb12 6 345.25 370.80 -166.62 333.25 49.346 1 2.146e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: largeDS2
## Models:
## modDSb6: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb6: (1 | Gsp)
## modDSb12: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb6 5 346.09 367.38 -168.04 336.09
## modDSb12 6 345.25 370.80 -166.62 333.25 2.8366 1 0.09214 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: largeDS2
## Models:
## modDSb4: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp)
## modDSb5: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSb5: (1 | Gsp)
## modDSb6: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb6: (1 | Gsp)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb4 3 350.54 363.31 -172.27 344.54
## modDSb5 4 345.35 362.38 -168.68 337.35 7.1909 1 0.007327 **
## modDSb6 5 346.09 367.38 -168.04 336.09 1.2629 1 0.261106
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: largeDS2
## Models:
## modDSb10: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDSb10: Ref)
## modDSb11: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSb11: (1 | Gsp) + (1 | Ref)
## modDSb12: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb10 4 352.54 369.57 -172.27 344.54
## modDSb11 5 343.95 365.24 -166.98 333.95 10.589 1 0.001137 **
## modDSb12 6 345.25 370.80 -166.62 333.25 0.701 1 0.402432
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: largeDS2
## Models:
## modDSb5: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSb5: (1 | Gsp)
## modDSb11: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSb11: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb5 4 345.35 362.38 -168.68 337.35
## modDSb11 5 343.95 365.24 -166.98 333.95 3.3985 1 0.06526 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Gsp) + (1 | Ref)
## Data: largeDS2
##
## AIC BIC logLik deviance df.resid
## 344.0 365.2 -167.0 334.0 517
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9660 -0.3198 -0.0383 -0.0265 4.5830
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 28.3671 5.3261
## Ref (Intercept) 0.9538 0.9766
## Number of obs: 522, groups: Gsp, 74; Ref, 18
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.702214 2.000071 -3.851 0.000118 ***
## Sed_level_stand_mg 0.003963 0.001377 2.878 0.003998 **
## Sed_exposure_days 0.026231 0.009663 2.715 0.006635 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_large_necroses ~ log10_conc + log10_dur + (1 | Gsp) +
## (1 | Ref)
## Data: largeDS2
##
## AIC BIC logLik deviance df.resid
## 322.9 344.2 -156.5 312.9 517
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7377 -0.1869 -0.0426 -0.0201 4.1335
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 27.088 5.205
## Ref (Intercept) 1.146 1.071
## Number of obs: 522, groups: Gsp, 74; Ref, 18
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -14.0241 2.7141 -5.167 2.38e-07 ***
## log10_conc 2.5721 0.5767 4.460 8.20e-06 ***
## log10_dur 2.7644 0.8983 3.077 0.00209 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Effect estimates
## [1] 0.9099617##
## p 0 1
## 0 434 34
## 1 13 41## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.9494## R2m R2c
## theoretical 0.09338887 0.9053851
## delta 0.09041040 0.8765094
## $Gsp
##
## $Ref
## Groups Name Std.Dev.
## Gsp (Intercept) 5.2046
## Ref (Intercept) 1.0707## [1] 2.917504
## Est LL UL
## (Intercept) 8.117388e-07 3.973205e-09 1.658409e-04
## log10_conc 1.309276e+01 4.227962e+00 4.054444e+01
## log10_dur 1.586911e+01 2.728155e+00 9.230725e+01For every 10-fold increase in exposure concentration of deposited sediment, the odds of any tissue mortality increase by 13.1 times (95% CI 4.2, 40.5; GLMM z = 4.460, p < 0.00001), while holding exposure duration constant and accounting for variability among studies and species. For every 10-fold change in exposure duration, the odds of any tissue mortality increase by 15.9 times (95% CI 2.7, 9.2; GLMM z = 3.077, p = 0.002), while holding exposure concentration constant and accounting for variability among studies and species.
Prediction figure
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDSb11_log <- predict(modDSb11_log, newdata=largeDS2, type="response")
largeDS3 <- cbind(largeDS2, pred_modDSb11_log)
largeDS_plot <- ggplot(data = largeDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_large_necroses,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Large necroses due to sediment exposure?",
color = "Study") +
scale_x_log10(limits=c(1,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDSb11_log), inherit.aes=FALSE) +
theme_bw()
largeDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
By exposure concentration
jvalues <- with(largeDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
largeDS2$log10_conc <- j
predict(modDSb11_log, newdata = largeDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
##
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
largeDS_plot2 <- ggplot() +
geom_point(data = largeDS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_large_necroses)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of large\nnecroses due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(1,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS2),
label=c("0.01","0.1","1","10","100","1000",round(loaelDS2,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelDS2, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
largeDS_plot2
By exposure duration
jvalues <- with(largeDS2, seq(from = min(log10_dur), to = max(log10_dur), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
largeDS2$log10_dur <- j
predict(modDSb11_log, newdata = largeDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
##
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_dur values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_dur", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
largeDS_plot2 <- ggplot() +
geom_point(data = largeDS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_large_necroses)) +
labs(x = expression("Sediment exposure duration (days)"),
y = "Predicted probability of large\nnecroses due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.1,max(largeDS2$Sed_exposure_days)), breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
largeDS_plot2
Total mortality
## n
## 1 509## Data: totalDS2
## Models:
## modDSc6: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSc6: (1 | Gsp)
## modDSc9: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSc9: (1 | Ref)
## modDSc12: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSc12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSc6 5 306.66 327.82 -148.33 296.66
## modDSc9 5 392.75 413.91 -191.37 382.75 0.000 0 1
## modDSc12 6 295.74 321.13 -141.87 283.74 99.009 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: totalDS2
## Models:
## modDSc10: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modDSc11: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSc11: (1 | Gsp) + (1 | Ref)
## modDSc12: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSc12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSc10 4 308.19 325.12 -150.09 300.19
## modDSc11 5 295.52 316.69 -142.76 285.52 14.664 1 0.0001285 ***
## modDSc12 6 295.74 321.13 -141.87 283.74 1.784 1 0.1816564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: totalDS2
## Models:
## modDSc4: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp)
## modDSc7: Binary_death_total ~ Sed_level_stand_mg + (1 | Ref)
## modDSc10: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSc4 3 319.21 331.90 -156.60 313.21
## modDSc7 3 396.32 409.02 -195.16 390.32 0.000 0 1
## modDSc10 4 308.19 325.12 -150.09 300.19 90.131 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_death_total ~ log10_conc + log10_dur + (1 | Gsp) + (1 |
## Ref)
## Data: totalDS2
##
## AIC BIC logLik deviance df.resid
## 255.0 276.2 -122.5 245.0 504
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2021 -0.0137 -0.0021 0.0000 7.6530
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 280.27 16.741
## Ref (Intercept) 42.47 6.517
## Number of obs: 509, groups: Gsp, 75; Ref, 23
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -38.909 6.755 -5.760 8.42e-09 ***
## log10_conc 6.103 1.500 4.069 4.73e-05 ***
## log10_dur 7.855 1.711 4.590 4.44e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Effect estimates
## [1] 0.9685658##
## p 0 1
## 0 419 9
## 1 7 74## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.9905## R2m R2c
## theoretical 0.0683368385 0.99059879
## delta 0.0009323897 0.01351576
## $Gsp
##
## $Ref
## Groups Name Std.Dev.
## Gsp (Intercept) 16.741
## Ref (Intercept) 6.517## [1] 676.5725
## Est LL UL
## (Intercept) 1.264562e-17 2.247393e-23 7.115428e-12
## log10_conc 4.471818e+02 2.363854e+01 8.459556e+03
## log10_dur 2.577456e+03 9.003714e+01 7.378375e+04For every 10-fold increase in exposure concentration of deposited sediment, the odds of total colony mortality increase by 447 times (95% CI 24, 8459; GLMM z = 4.069, p < 0.0001), after accounting for exposure duration and variability among studies and species. For every 10-fold change in exposure duration, the odds of total colony mortality increase by 2577 times (95% CI 90, 73,784; GLMM z = 4.590, p = 0.00001), after accounting for exposure concentration and variability among studies and species.
Prediction figure
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDSc11_log <- predict(modDSc11_log, newdata=totalDS2, type="response")
totalDS3 <- cbind(totalDS2, pred_modDSc11_log)
totalDS_plot <- ggplot(data = totalDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_death_total,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Total mortality due to sediment exposure?",
color = "Study") +
scale_x_log10(limits=c(0.01,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDSc11_log), inherit.aes=FALSE) +
theme_bw()
totalDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
By exposure concentration
jvalues <- with(totalDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
totalDS2$log10_conc <- j
predict(modDSc11_log, newdata = totalDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
##
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
totalDS_plot2 <- ggplot() +
geom_point(data = totalDS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_death_total,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of total colony\nmortality due to sediment exposure",
color = "Study",
linetype = "Predicted Probability") +
scale_x_log10(limits=c(0.01,1000), breaks=c(0.01,0.1,1,10,100,1000,loaelDS3),
label=c("0.01","0.1","1","10","100","1000",round(loaelDS3,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelDS3, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
totalDS_plot2
By exposure duration
jvalues <- with(totalDS2, seq(from = min(log10_dur), to = max(log10_dur), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
totalDS2$log10_dur <- j
predict(modDSc11_log, newdata = totalDS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_dur values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_dur", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
totalDS_plot2 <- ggplot() +
geom_point(data = totalDS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_death_total)) +
labs(x = expression("Sediment exposure duration (days)"),
y = "Predicted probability of total colony\nmortality due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.1,max(totalDS2$Sed_exposure_days)), breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
totalDS_plot2
SUSPENDED SEDIMENT
Small Necroses
## n
## 1 147## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Data: smallSS2
## Models:
## modSS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## modSS9: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS9: (1 | Ref)
## modSS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## modSS15: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS15: (1 | Ref_name/Ref)
## modSS18: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## modSS21: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS21: (1 | Gsp) + (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS6 5 108.19 123.14 -49.094 98.189
## modSS9 5 126.29 141.24 -58.145 116.290 0.000 0 1
## modSS12 6 110.19 128.13 -49.094 98.189 18.102 1 2.094e-05 ***
## modSS15 6 128.29 146.23 -58.145 116.290 0.000 0 1
## modSS18 6 110.19 128.13 -49.094 98.189 18.102 0 < 2.2e-16 ***
## modSS21 7 112.19 133.12 -49.094 98.189 0.000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: smallSS2
## Models:
## modSS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## modSS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## modSS18: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS6 5 108.19 123.14 -49.094 98.189
## modSS12 6 110.19 128.13 -49.094 98.189 0 1 1
## modSS18 6 110.19 128.13 -49.094 98.189 0 0 1## Data: smallSS2
## Models:
## modSS4: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp)
## modSS5: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS5: (1 | Gsp)
## modSS6: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS4 3 135.52 144.49 -64.759 129.518
## modSS5 4 106.75 118.71 -49.375 98.750 30.7680 1 2.908e-08 ***
## modSS6 5 108.19 123.14 -49.094 98.189 0.5607 1 0.454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: smallSS2
## Models:
## modSS10: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSS10: Ref)
## modSS11: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS11: (1 | Gsp) + (1 | Ref)
## modSS12: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS10 4 137.45 149.41 -64.724 129.447
## modSS11 5 108.75 123.70 -49.375 98.750 30.6976 1 3.015e-08 ***
## modSS12 6 110.19 128.13 -49.094 98.189 0.5607 1 0.454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: smallSS2
## Models:
## modSS16: Binary_small_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSS16: Ref_name)
## modSS17: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS17: (1 | Gsp) + (1 | Ref_name)
## modSS18: Binary_small_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS16 4 137.45 149.41 -64.724 129.447
## modSS17 5 108.75 123.70 -49.375 98.750 30.6976 1 3.015e-08 ***
## modSS18 6 110.19 128.13 -49.094 98.189 0.5607 1 0.454
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: smallSS2
## Models:
## modSS5: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS5: (1 | Gsp)
## modSS11: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS11: (1 | Gsp) + (1 | Ref)
## modSS17: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS17: (1 | Gsp) + (1 | Ref_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS5 4 106.75 118.71 -49.375 98.75
## modSS11 5 108.75 123.70 -49.375 98.75 0 1 1
## modSS17 5 108.75 123.70 -49.375 98.75 0 0 1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_small_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Gsp)
## Data: smallSS2
##
## AIC BIC logLik deviance df.resid
## 106.7 118.7 -49.4 98.7 143
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9893 -0.3421 -0.1326 -0.0291 3.8738
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 3.44 1.855
## Number of obs: 147, groups: Gsp, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.06121 1.75339 -4.598 4.28e-06 ***
## Sed_level_stand_mg 0.03325 0.00998 3.331 0.000864 ***
## Sed_exposure_days 0.07866 0.01806 4.355 1.33e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Effect estimates
## [1] 0.8707483##
## p 0 1
## 0 108 11
## 1 8 20## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.9189## R2m R2c
## theoretical 0.3536090 0.6839914
## delta 0.2701203 0.5224978
## $Gsp
## Groups Name Std.Dev.
## Gsp (Intercept) 1.8546## [1] 6.389086
## Est LL UL
## (Intercept) 0.003743975 0.0003457821 0.04053811
## Sed_level_stand_mg 1.023311817 1.0095310571 1.03728069
## Sed_exposure_days 1.056033682 1.0304395187 1.08226356## Est LL UL
## (Intercept) 8.685337e-09 3.177978e-12 2.373682e-05
## Sed_level_stand_mg 1.079558e+00 1.032013e+00 1.129292e+00
## Sed_exposure_days 1.198549e+00 1.104739e+00 1.300325e+00For every 10-fold increase in exposure concentration of suspended sediment, the odds of developing small tissue necroses increase by 1.08 times (95% CI 1.03, 1.13; GLMM z = 3.331, p = 0.0008), while holding exposure duration constant and after accounting for variability among species. For every 10-fold increase in exposure duration, the odds of small tissue necroses increase by 1.20 times (95% CI 1.10, 1.30; GLMM z = 4.355, p < 0.0001), while holding exposure concentration constant and after accounting for variability among species.
Prediction figure
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSS5 <- predict(modSS5, newdata=smallSS2, type="response")
smallSS3 <- cbind(smallSS2, pred_modSS5)
smallSS_plot <- ggplot(data = smallSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_small_necroses,
color = Gsp)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Small necroses due to sediment exposure?",
color = "Species") +
scale_x_log10(limits=c(1,max(smallSS2$Sed_level_stand_mg)),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSS5), inherit.aes=FALSE) +
theme_bw()
smallSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
By sediment exposure concentration
jvalues <- with(smallSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
smallSS2$Sed_level_stand_mg <- j
predict(modSS5, newdata = smallSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
smallSS_plot2 <- ggplot() +
geom_point(data = smallSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_small_necroses)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of small\nnecroses due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,125), breaks=c(0,25,50,75,100,125,loaelSS1),
labels=c("0","25","50","75","100","125",round(loaelSS1,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS1, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
smallSS_plot2
By sediment exposure duration
jvalues <- with(smallSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
smallSS2$Sed_exposure_days <- j
predict(modSS5, newdata = smallSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
smallSS_plot3 <- ggplot() +
geom_point(data = smallSS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_small_necroses,
color = Gsp)) +
labs(x = "Sediment exposure duration (days)",
y = "Predicted probability of small\nnecroses to sediment exposure",
color = "Species",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,max(smallSS2$Sed_exposure_days))) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
smallSS_plot3
Large Necroses
## n
## 1 147## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## df AIC
## modSSb1 2 52.19214
## modSSb2 3 11.54518
## modSSb3 4 13.54518## Data: largeSS2
## Models:
## modSSb7: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Ref)
## modSSb8: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb8: (1 | Ref)
## modSSb9: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSb9: (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb7 3 31.918 40.889 -12.9588 25.9176
## modSSb8 4 13.545 25.507 -2.7726 5.5452 20.372 1 6.374e-06 ***
## modSSb9 5 15.545 30.497 -2.7726 5.5452 0.000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: largeSS2
## Models:
## modSSb13: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modSSb14: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb14: (1 | Ref_name/Ref)
## modSSb15: Binary_large_necroses ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSb15: (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb13 4 33.918 45.879 -12.9588 25.9176
## modSSb14 5 15.545 30.497 -2.7726 5.5452 20.372 1 6.374e-06 ***
## modSSb15 6 17.545 35.488 -2.7726 5.5452 0.000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: largeSS2
## Models:
## modSSb4: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp)
## modSSb8: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb8: (1 | Ref)
## modSSb10: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSSb10: Ref)
## modSSb16: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSSb16: Ref_name)
## modSSb14: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb14: (1 | Ref_name/Ref)
## modSSb19: Binary_large_necroses ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSSb19: Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb4 3 36.784 45.756 -15.3921 30.7842
## modSSb8 4 13.545 25.507 -2.7726 5.5452 25.239 1 5.065e-07 ***
## modSSb10 4 33.803 45.765 -12.9014 25.8028 0.000 0 1
## modSSb16 4 33.803 45.765 -12.9014 25.8028 0.000 0 1
## modSSb14 5 15.545 30.497 -2.7726 5.5452 20.258 1 6.768e-06 ***
## modSSb19 5 35.803 50.755 -12.9014 25.8028 0.000 0 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: largeSS2
## Models:
## modSSb8: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb8: (1 | Ref)
## modSSb14: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSb14: (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb8 4 13.545 25.507 -2.7726 5.5452
## modSSb14 5 15.545 30.497 -2.7726 5.5452 0 1 1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Ref)
## Data: largeSS2
##
## AIC BIC logLik deviance df.resid
## 13.5 25.5 -2.8 5.5 143
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1 0 0 0 1
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 0 0
## Number of obs: 147, groups: Ref, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -779.798 968.212 -0.805 0.4206
## Sed_level_stand_mg 2.026 31.654 0.064 0.9490
## Sed_exposure_days 8.582 3.984 2.154 0.0312 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular## boundary (singular) fit: see ?isSingular## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_large_necroses ~ log10_conc + log10_dur + (1 | Ref)
## Data: largeSS2
##
## AIC BIC logLik deviance df.resid
## 13.5 25.5 -2.8 5.5 143
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1 0 0 0 1
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 3.193e-17 5.65e-09
## Number of obs: 147, groups: Ref, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1347.0 49256351.2 0 1
## log10_conc 85.8 16165394.5 0 1
## log10_dur 634.7 15817710.9 0 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular##
## Call:
## glm(formula = Binary_large_necroses ~ Sed_level_stand_mg + Sed_exposure_days,
## family = binomial(link = "logit"), data = largeSS2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.177 0.000 0.000 0.000 1.177
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -390.084 100469.392 -0.004 0.997
## Sed_level_stand_mg 1.077 314.538 0.003 0.997
## Sed_exposure_days 4.271 1100.945 0.004 0.997
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 50.1358 on 146 degrees of freedom
## Residual deviance: 5.5452 on 144 degrees of freedom
## AIC: 11.545
##
## Number of Fisher Scoring iterations: 25Effect estimates
## [1] 0.9863946##
## p 0 1
## 0 141 2
## 1 0 4## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.9976## R2m R2c
## theoretical 0.9996563 0.9996563
## delta 0.9973381 0.9973381There is no significant relationship between sediment exposure and the odds of developing large tissue necroses (concentration: GLMM z = 0.003, p = 0.997; duration: GLMM z = 0.004, p = 0.997).
Prediction figure
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSb2 <- predict(modSSb2, newdata=largeSS2, type="response")
largeSS3 <- cbind(largeSS2, pred_modSSb2)
largeSS_plot <- ggplot(data = largeSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_large_necroses,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Large necroses due to sediment exposure?",
color = "Study") +
scale_x_continuous(limits=c(0,max(largeSS2$Sed_level_stand_mg))) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSb2), inherit.aes=FALSE) +
theme_bw()
largeSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
By exposure concentration
jvalues <- with(largeSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
largeSS2$Sed_level_stand_mg <- j
predict(modSSb2, newdata = largeSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
largeSS_plot2 <- ggplot() +
geom_point(data = largeSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_large_necroses,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of large\nnecroses due to sediment exposure",
color = "Binary data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,max(largeSS2$Sed_level_stand_mg)),
breaks = c(0,25,50,75,100,125,loaelSS2),
labels = c("0","25","50","75","100","125",round(loaelSS2,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS2, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
largeSS_plot2
By exposure duration
jvalues <- with(largeSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
largeSS2$Sed_exposure_days <- j
predict(modSSb2, newdata = largeSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
largeSS_plot2 <- ggplot() +
geom_point(data = largeSS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_large_necroses)) +
labs(x = expression("Sediment exposure duration (days)"),
y = "Predicted probability of large\nnecroses due to sediment exposure",
linetype = "Predicted\nProbability") +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
largeSS_plot2
Total mortality
## n
## 1 176## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## df AIC
## modSSc1 2 56.11582
## modSSc2 3 11.54518
## modSSc3 4 13.54518## Data: totalSS2
## Models:
## modSSc7: Binary_death_total ~ Sed_level_stand_mg + (1 | Ref)
## modSSc8: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSc8: (1 | Ref)
## modSSc9: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSc9: (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSc7 3 32.442 41.953 -13.2208 26.4416
## modSSc8 4 13.545 26.227 -2.7726 5.5452 20.896 1 4.848e-06 ***
## modSSc9 5 15.545 31.398 -2.7726 5.5452 0.000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: totalSS2
## Models:
## modSSc13: Binary_death_total ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modSSc15: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSc15: (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSc13 4 34.441 47.122 -13.2203 26.4406
## modSSc15 6 17.545 36.568 -2.7726 5.5452 20.895 2 2.902e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: totalSS2
## Models:
## modSSc4: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp)
## modSSc8: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSc8: (1 | Ref)
## modSSc10: Binary_death_total ~ Sed_level_stand_mg + (1 | Gsp) + (1 | Ref)
## modSSc15: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSc15: (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSc4 3 38.093 47.605 -16.0466 32.093
## modSSc8 4 13.545 26.227 -2.7726 5.545 26.548 1 2.571e-07 ***
## modSSc10 4 34.442 47.124 -13.2208 26.442 0.000 0 1
## modSSc15 6 17.545 36.568 -2.7726 5.545 20.896 2 2.900e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: totalSS2
## Models:
## modSSc8: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +
## modSSc8: (1 | Ref)
## modSSc15: Binary_death_total ~ Sed_level_stand_mg * Sed_exposure_days +
## modSSc15: (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSc8 4 13.545 26.227 -2.7726 5.5452
## modSSc15 6 17.545 36.568 -2.7726 5.5452 0 2 1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Ref)
## Data: totalSS2
##
## AIC BIC logLik deviance df.resid
## 13.5 26.2 -2.8 5.5 172
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1 0 0 0 1
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 5.443e-17 7.378e-09
## Number of obs: 176, groups: Ref, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -862.270 967.393 -0.891 0.3728
## Sed_level_stand_mg 2.029 31.628 0.064 0.9489
## Sed_exposure_days 9.562 3.984 2.400 0.0164 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## boundary (singular) fit: see ?isSingular## boundary (singular) fit: see ?isSingular## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_death_total ~ log10_conc + log10_dur + (1 | Ref)
## Data: totalSS2
##
## AIC BIC logLik deviance df.resid
## 13.5 26.2 -2.8 5.5 172
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1 0 0 0 1
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 0 0
## Number of obs: 176, groups: Ref, 8
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2676.5 2071.5 -1.292 0.196
## log10_conc 85.2 4512.3 0.019 0.985
## log10_dur 1326.1 3518.0 0.377 0.706
## convergence code: 0
## boundary (singular) fit: see ?isSingular##
## Call:
## glm(formula = Binary_death_total ~ Sed_level_stand_mg + Sed_exposure_days,
## family = binomial(link = "logit"), data = totalSS2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.177 0.000 0.000 0.000 1.177
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -385.217 88625.663 -0.004 0.997
## Sed_level_stand_mg 1.064 281.353 0.004 0.997
## Sed_exposure_days 4.217 968.684 0.004 0.997
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 52.3378 on 175 degrees of freedom
## Residual deviance: 5.5452 on 173 degrees of freedom
## AIC: 11.545
##
## Number of Fisher Scoring iterations: 25Effect estimates
## [1] 0.9886364##
## p 0 1
## 0 170 2
## 1 0 4## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.998## R2m R2c
## theoretical 0.9996362 0.9996362
## delta 0.9966523 0.9966523There is no significant relationship between sediment exposure and the odds of total colony mortality (concentration: GLMM z = 0.004, p = 0.997; duration: GLMM z = 0.004, p = 0.997).
Prediction figure
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSc2 <- predict(modSSc2, newdata=totalSS2, type="response")
totalSS3 <- cbind(totalSS2, pred_modSSc2)
totalSS_plot <- ggplot(data = totalSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_death_total,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Total mortality due to sediment exposure?",
color = "Study") +
scale_x_log10(limits=c(1,200),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSc2), inherit.aes=FALSE) +
theme_bw()
totalSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
By exposure concentration
jvalues <- with(totalSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
totalSS2$Sed_level_stand_mg <- j
predict(modSSc2, newdata = totalSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
totalSS_plot2 <- ggplot() +
geom_point(data = totalSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_death_total,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of total colony\nmortality due to sediment exposure",
color = "Binary data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(1,200), breaks=c(0.01,0.1,1,10,100,1000,loaelSS3),
label=c("0.01","0.1","1","10","100","1000",round(loaelSS3,digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS3, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
totalSS_plot2
By exposure duration
jvalues <- with(totalSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
totalSS2$Sed_exposure_days <- j
predict(modSSc2, newdata = totalSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
totalSS_plot2 <- ggplot() +
geom_point(data = totalSS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_death_total,
color = Ref)) +
labs(x = expression("Sediment exposure duration (days)"),
y = "Predicted probability of total colony\nmortality due to sediment exposure",
color = "Study",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,max(totalSS2$Sed_exposure_days))) +
geom_ribbon(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
totalSS_plot2